clear
set more off
cd "SET YOUR DIRECTORY HERE"
set seed 6451166

	
************************
***Table 1 & Figure 2***		
************************

clear
use "Mass survey.dta"
set more off

	*set up variables that we will use during the loops
		gen pred1 = .
		gen pred2 = .
		gen pred3 = .
		gen pred4 = .
		gen pred5 = .

		capture drop xvalue
		gen xvalue = . 

		local a = 1	
		
		while `a' <= 7 {
			replace xvalue = `a' in `a'
			
			local a = `a' + 1						
		}				
	
	
	*roads and paving
		meoprobit roads_and_paving tamciud income gvtidentification cdummy* 	
		meoprobit roads_and_paving tamciud cdummy*
	
		local a = 1
			
		while `a' <= 7 {
			
			replace pred1 = normprob(_b[/cut1] - (_b[tamciud]*`a')) in `a'
			replace pred2 = (normprob(_b[/cut2] - (_b[tamciud]*`a'))) - (normprob(_b[/cut1] - (_b[tamciud]*`a'))) in `a'
			replace pred3 = (normprob(_b[/cut3] - (_b[tamciud]*`a'))) - (normprob(_b[/cut2] - (_b[tamciud]*`a'))) in `a'
			replace pred4 = (normprob(_b[/cut4] - (_b[tamciud]*`a'))) - (normprob(_b[/cut3] - (_b[tamciud]*`a'))) in `a'
			replace pred5 = 1-(normprob(_b[/cut4] - (_b[tamciud]*`a'))) in `a'	
								
			local a = `a' + 1					/* generate the next value for the xaxis */			
				
			}		
					
		*Get cumultative predicted probabilities 	
			replace pred2 = pred1 + pred2
			replace pred3 = pred2 + pred3
			replace pred4 = pred3 + pred4
			replace pred5 = pred4 + pred5
		
		twoway (area pred1 xvalue,  lcolor(black)) (rarea pred1 pred2 xvalue,  lcolor(black)) (rarea pred2 pred3 xvalue,  lcolor(black)) (rarea pred3 pred4 xvalue,  lcolor(black)) (rarea pred4 pred5 xvalue,  lcolor(black)), legend(order(1 "Very satisfied" 2 "Rather satisfied" 3 "Not very satisfied" 4 "Not at all satisfied" 5 "Doesn't have these services")) name(tempg1, replace) graphregion(color(white)) title("Roads and Paving") ytitle("") xlabel(1(1)7) ylab(0(0.2)1, nogrid) xtitle("") ytitle("Cumulative probability") xlabel(1 "Small" 2 "2" 3 "3" 4 "4" 5 "5" 6 "6" 7 "Large")

	*Public transportation
		meoprobit public_transport tamciud income gvtidentification cdummy* 	
		meoprobit public_transport tamciud cdummy*
	
		local a = 1
			
		while `a' <= 7 {
			
			replace pred1 = normprob(_b[/cut1] - (_b[tamciud]*`a')) in `a'
			replace pred2 = (normprob(_b[/cut2] - (_b[tamciud]*`a'))) - (normprob(_b[/cut1] - (_b[tamciud]*`a'))) in `a'
			replace pred3 = (normprob(_b[/cut3] - (_b[tamciud]*`a'))) - (normprob(_b[/cut2] - (_b[tamciud]*`a'))) in `a'
			replace pred4 = (normprob(_b[/cut4] - (_b[tamciud]*`a'))) - (normprob(_b[/cut3] - (_b[tamciud]*`a'))) in `a'
			replace pred5 = 1-(normprob(_b[/cut4] - (_b[tamciud]*`a'))) in `a'	
								
			local a = `a' + 1					/* generate the next value for the xaxis */			
				
			}		
				
		*Get cumultative predicted probabilities 	
			replace pred2 = pred1 + pred2
			replace pred3 = pred2 + pred3
			replace pred4 = pred3 + pred4
			replace pred5 = pred4 + pred5
		
		twoway (area pred1 xvalue,  lcolor(black)) (rarea pred1 pred2 xvalue,  lcolor(black)) (rarea pred2 pred3 xvalue,  lcolor(black)) (rarea pred3 pred4 xvalue,  lcolor(black)) (rarea pred4 pred5 xvalue,  lcolor(black)), legend(order(1 "Very satisfied" 2 "Rather satisfied" 3 "Not very satisfied" 4 "Not at all satisfied" 5 "Doesn't have these services")) name(tempg2, replace) graphregion(color(white)) title("Public Transportation") ytitle("") xlabel(1(1)7) ylab(0(0.2)1, nogrid) xtitle("") xlabel(1 "Small" 2 "2" 3 "3" 4 "4" 5 "5" 6 "6" 7 "Large")
		
	*Garbage
		meoprobit garbage tamciud income gvtidentification cdummy* 	
		meoprobit garbage tamciud cdummy*
	
		local a = 1
			
		while `a' <= 7 {
			
			replace pred1 = normprob(_b[/cut1] - (_b[tamciud]*`a')) in `a'
			replace pred2 = (normprob(_b[/cut2] - (_b[tamciud]*`a'))) - (normprob(_b[/cut1] - (_b[tamciud]*`a'))) in `a'
			replace pred3 = (normprob(_b[/cut3] - (_b[tamciud]*`a'))) - (normprob(_b[/cut2] - (_b[tamciud]*`a'))) in `a'
			replace pred4 = (normprob(_b[/cut4] - (_b[tamciud]*`a'))) - (normprob(_b[/cut3] - (_b[tamciud]*`a'))) in `a'
			replace pred5 = 1-(normprob(_b[/cut4] - (_b[tamciud]*`a'))) in `a'	
								
			local a = `a' + 1					/* generate the next value for the xaxis */			
				
			}		
				
				
		*Get cumultative predicted probabilities 	
			replace pred2 = pred1 + pred2
			replace pred3 = pred2 + pred3
			replace pred4 = pred3 + pred4
			replace pred5 = pred4 + pred5
		
		twoway (area pred1 xvalue,  lcolor(black)) (rarea pred1 pred2 xvalue,  lcolor(black)) (rarea pred2 pred3 xvalue,  lcolor(black)) (rarea pred3 pred4 xvalue,  lcolor(black)) (rarea pred4 pred5 xvalue,  lcolor(black)), legend(order(1 "Very satisfied" 2 "Rather satisfied" 3 "Not very satisfied" 4 "Not at all satisfied" 5 "Doesn't have these services")) name(tempg3, replace) graphregion(color(white)) title("Garbage") ytitle("") xlabel(1(1)7) ylab(0(0.2)1, nogrid) ytitle("Cumulative probability") xtitle("Size of town/city") xlabel(1 "Small" 2 "2" 3 "3" 4 "4" 5 "5" 6 "6" 7 "Large")
		
	*Sewerage
		meoprobit sewerage tamciud income gvtidentification cdummy* 	
		meoprobit sewerage tamciud cdummy*
	
		local a = 1
			
		while `a' <= 7 {
			
			replace pred1 = normprob(_b[/cut1] - (_b[tamciud]*`a')) in `a'
			replace pred2 = (normprob(_b[/cut2] - (_b[tamciud]*`a'))) - (normprob(_b[/cut1] - (_b[tamciud]*`a'))) in `a'
			replace pred3 = (normprob(_b[/cut3] - (_b[tamciud]*`a'))) - (normprob(_b[/cut2] - (_b[tamciud]*`a'))) in `a'
			replace pred4 = (normprob(_b[/cut4] - (_b[tamciud]*`a'))) - (normprob(_b[/cut3] - (_b[tamciud]*`a'))) in `a'
			replace pred5 = 1-(normprob(_b[/cut4] - (_b[tamciud]*`a'))) in `a'	
								
			local a = `a' + 1					/* generate the next value for the xaxis */			
				
			}		
				
		*Get cumultative predicted probabilities 	
			replace pred2 = pred1 + pred2
			replace pred3 = pred2 + pred3
			replace pred4 = pred3 + pred4
			replace pred5 = pred4 + pred5
		
		twoway (area pred1 xvalue,  lcolor(black)) (rarea pred1 pred2 xvalue,  lcolor(black)) (rarea pred2 pred3 xvalue,  lcolor(black)) (rarea pred3 pred4 xvalue,  lcolor(black)) (rarea pred4 pred5 xvalue,  lcolor(black)), legend(order(1 "Very satisfied" 2 "Rather satisfied" 3 "Not very satisfied" 4 "Not at all satisfied" 5 "Doesn't have these services")) name(tempg4, replace) graphregion(color(white)) title("Sewerage") ytitle("") xlabel(1(1)7) ylab(0(0.2)1, nogrid) xtitle("Size of town/city") xlabel(1 "Small" 2 "2" 3 "3" 4 "4" 5 "5" 6 "" 7 "Large")
	
		grc1leg tempg1 tempg2 tempg3 tempg4, graphregion(color(white))

************************
***Table 2 & Figure 3***
************************

clear
use "Mass survey.dta"
set more off
	
	*Deputy contact
		meoprobit deputy_contact tamciud  cdummy*
		meoprobit deputy_contact tamciud income gvtidentification cdummy* 	
		
	*Demand	
		probit index tamciud income gvtidentification cdummy*
		probit index tamciud cdummy*	

		
		do "simUtils"

		  tempname b V sig dfsig
	      matrix `b' = e(b)                            /* 1 x k vector         */
	      matrix `V' = e(V)                            /* k x k variance matrix*/
	      local N = `e(N)'                             /* save # observations  */
	
	      _simp, b(`b') v(`V') s(10000) g(bbt)		/* run the program to generate the simulations */ 
	

			generate predVAL=.			//Variables that I will fill in during the loop
			generate predHIGH=.			//Variables that I will fill in during the loop
			generate predLOW=.			//Variables that I will fill in during the loop
			generate xaxis=.			//Variables that I will fill in during the loop

		
			local a=1			
			local b=1			
		
			while `a' <= 7 {
				replace xaxis = `a' in `b'
		
				local MRP = `a'
				
				capture drop pred
				gen pred = normprob(bbt15 + bbt1*`a')
	
				_pctile pred, p(2.5,97.5)  			/* stores the 5 and 95 percentiles of pred under the 
														internal variables r(r1) and r(r2) */

				replace predLOW = r(r1) in `b'		/* assigns the value of the 5 percetile to the variable predLOW */
		
				replace predHIGH = r(r2) in `b' 	/* assigns the value of the 95 percetile to the variable predHIGH */

				summarize pred						/* stores the mean of pred in internal variable r(mean) */
				replace predVAL = r(mean) in `b'
				capture drop pred
				
				local a = `a' + 1					/* generate the next value for the xaxis */
				local b = `b' + 1					/* iterate the loop */

	}							
		
		twoway (line predVAL xaxis, lcolor(black)) (line predLOW xaxis, lcolor(black) lpattern(dash)) (line predHIGH xaxis, lcolor(black) lpattern(dash)), xtitle("")  ytitle("Pr(Relative Deprivation and Contact)") legend(off) graphregion(color(white)) ylabel(,nogrid) xtitle("Size of town/city")  xlabel(1 "Small" 2 "2" 3 "3" 4 "4" 5 "5" 6 "6" 7 "Large") 	

************************
***Table 3 & Figure 4***		
************************
clear
use "Mass survey.dta"
set more off

	*merge data to get population density and other variables
		merge m:m district_name session_beg using "Elite survey"
		keep if _merge == 3
	
	*drop districts with only very few observations
		egen districtobs = count(numinves), by (district_name)
		drop if districtobs < 40
		
	*run analysis
		collapse index idenpa perc_vote gdp_pc_region2 district_pop_density2 cdummy*, by(district_name)			
		reg index district_pop_density2 gdp_pc_region2 perc_vote  cdummy*
		reg index district_pop_density2 cdummy* 
	
	*Figure
		do "simUtils"

		  tempname b V sig dfsig
	      matrix `b' = e(b)                            /* 1 x k vector         */
	      matrix `V' = e(V)                            /* k x k variance matrix*/
	      local N = `e(N)'                             /* save # observations  */
	
	      _simp, b(`b') v(`V') s(10000) g(bbt)		/* run the program to generate the simulations */ 
	

			generate predVAL=.			//variable that I will fill in during the loop
			generate predHIGH=.			//variable that I will fill in during the loop
			generate predLOW=.			//variable that I will fill in during the loop
			generate xaxis=.			//variable that I will fill in during the loop

		
			local a=2.45824				//10th percentile
			local b=1			
		
			while `a' <= 5.837663 {
				replace xaxis = `a' in `b'
		
				local MRP = `a'
				
				capture drop pred
				gen pred = (bbt15 + (bbt2 + bbt5 + bbt6 + bbt7 + bbt8 + bbt9 + bbt10 + bbt11 + bbt12 + bbt13)/10) + bbt1*`a'
	
				_pctile pred, p(2.5,97.5)  			/* stores the 5 and 95 percentiles of pred under the 
														internal variables r(r1) and r(r2) */

				replace predLOW = r(r1) in `b'		/* assigns the value of the 5 percetile to the variable predLOW */
		
				replace predHIGH = r(r2) in `b' 	/* assigns the value of the 95 percetile to the variable predHIGH */

				summarize pred						/* stores the mean of pred in internal variable r(mean) */
				replace predVAL = r(mean) in `b'
				capture drop pred
				
				local a = `a' + 0.01				/* generate the next value for the xaxis */
				local b = `b' + 1					/* iterate the loop */

	}							
		replace predLOW = . if predLOW < 0
		replace xaxis = exp(xaxis)				//delog population density
		twoway (line predVAL xaxis, xscale(range(0 350)) xlabel(0(100)300) lcolor(black)) (line predLOW xaxis, lcolor(black) lpattern(dash)) (line predHIGH xaxis, lcolor(black) lpattern(dash)), xtitle("Population density in square kilometers")  ytitle("Percentage Demand") legend(off) graphregion(color(white)) ylabel(,nogrid)  

************************
***Table 4 & Figure 5***		
************************

clear
use "Elite survey.dta"
set more off

	*run analysis with fixed effects
		egen groupindicator = group(ccode survey_wave)
		egen districtindicator = group(ccode district_name)
		tab groupindicator, g(groupdummy)		
		
		*run model
			meoprobit q_24_01 district_pop_density2 dm2 government gdp_pc_region2 groupdummy*  || districtindicator:	
			meoprobit q_24_01 district_pop_density2 dm2 government groupdummy*  || districtindicator:	
				
		*create a pretty figure 
			gen pred1 = .
			gen pred2 = .
			gen pred3 = .
			gen pred4 = .
			gen pred5 = .

			capture drop xvalue
			gen xvalue = . 

				local a = -2
				local b = 1
				
				while `a' <= 10 {
				
				replace pred1 = normprob(_b[/cut1] - (_b[district_pop_density2]*`a' + _b[dm2]*1.791759)) in `b'
				replace pred2 = (normprob(_b[/cut2] - (_b[district_pop_density2]*`a' + _b[dm2]*1.791759))) - (normprob(_b[/cut1] - (_b[district_pop_density2]*`a' + _b[dm2]*1.791759))) in `b'
				replace pred3 = (normprob(_b[/cut3] - (_b[district_pop_density2]*`a' + _b[dm2]*1.791759))) - (normprob(_b[/cut2] - (_b[district_pop_density2]*`a' + _b[dm2]*1.791759))) in `b'
				replace pred4 = 1-(normprob(_b[/cut3] - (_b[district_pop_density2]*`a' + _b[dm2]*1.791759))) in `b'		
				replace xvalue = `a' in `b'
				local a = `a' + 1				
				local b = `b' + 1			
				}		
					
					
			*Get cumultative predicted probabilities 	
				replace pred2 = pred1 + pred2
				replace pred3 = pred2 + pred3
				replace pred4 = pred3 + pred4
			
			*gen figure
				twoway (area pred1 xvalue,  lcolor(black)) (rarea pred1 pred2 xvalue,  lcolor(black)) (rarea pred2 pred3 xvalue,  lcolor(black)) (rarea pred3 pred4 xvalue,  lcolor(black)), legend(order(1 "Not at all important" 2 "Little importance" 3 "Fairly important" 4 "Very important")) name(tempg1, replace) graphregion(color(white)) ytitle("Cumulative probability") xlabel(-2(2)10) ylab(0(0.2)1, nogrid) xtitle("Ln(Population density)") 
				
************************
***Table 5 & Figure 6***		
************************		

clear
use "Mass survey.dta"
set more off
set seed 8329793	

	*drop districts that were very small
		egen districtobs = count(numinves), by (district_name)
		drop if districtobs < 40
		
	*generate percentage of respondents who had an index value of "1" in each district
		egen index_dist = mean(index), by(district_name)
		
		keep session_beg index_dist district_name
		collapse session_beg index_dist, by(district_name)

	*save mass survey data
		saveold "indexdata.dta", replace

	*get elite data
		clear
		use "Elite survey.dta"
		set more off
		
	*merge to get index variable	
		merge m:1 district_name using "indexdata"
	
	*keep only countries that we have data for around the time of 2008
		keep if cou_name == "ARGENTINA" | cou_name == "COSTA RICA" | cou_name == "DOMINICAN REPUBLIC" | cou_name == "ECUADOR" | cou_name == "EL SALVADOR" | cou_name == "GUATEMALA" | cou_name == "HONDURAS" | cou_name == "NICARAGUA" | cou_name == "PARAGUAY" | cou_name == "PERU" | cou_name == "URUGUAY" 
		keep if survey_year > 2002 & survey_year < 2014 & survey_year ~= . 	
		
	*regression, first stage
		tab cou_name, g(cdummy)

		capture program drop firstststage
		program firstststage
			reg index_dist district_pop_density2 dm2 government cdummy1 cdummy2 cdummy3 cdummy4 cdummy5 cdummy6 cdummy7 cdummy8 cdummy9 cdummy10  if q_24_01 ~= . 
		end
		bootstrap, reps(100): firstststage
		
	*2. stage
		capture program drop my2sls
		program my2sls
		reg index_dist district_pop_density2 dm2 government cdummy1 cdummy2 cdummy3 cdummy4 cdummy5 cdummy6 cdummy7 cdummy8 cdummy9 cdummy10  if q_24_01 ~= . 
			predict IV if e(sample), xb
			oprobit q_24_01 IV dm2 government cdummy1 cdummy2 cdummy3 cdummy4 cdummy5 cdummy6 cdummy7 cdummy8 cdummy9 cdummy10 
			drop IV
		end
		bootstrap, reps(100): my2sls		

			
	*show figure 
		reg index_dist district_pop_density2 dm2 government cdummy1 cdummy2 cdummy3 cdummy4 cdummy5 cdummy6 cdummy7 cdummy8 cdummy9 cdummy10  if q_24_01 ~= . 
		predict IV if e(sample), xb
		oprobit q_24_01 IV dm2 government cdummy1 cdummy2 cdummy3 cdummy4 cdummy5 cdummy6 cdummy7 cdummy8 cdummy9 cdummy10  

		gen pred1 = .
		gen pred2 = .
		gen pred3 = .
		gen pred4 = .
		gen pred5 = .

		capture drop xvalue
		gen xvalue = . 

		local a = 0	
		local b = 1			
		
		while `a' <= 59 {
			replace xvalue = `a' in `b'
			
			local a = `a' + 1		
			local b = `b' + 1		
		}				
		
				local a = 0
				local b = 1
				
				while `a' <= 59 {
				
			replace pred1 = normprob(_b[/cut1] - (_b[IV]*`a' + _b[dm2]*2.302585)) in `b'
				replace pred2 = (normprob(_b[/cut2] - (_b[IV]*`a' + _b[dm2]*2.302585))) - (normprob(_b[/cut1] - (_b[IV]*`a' + _b[dm2]*2.302585))) in `b'
				replace pred3 = (normprob(_b[/cut3] - (_b[IV]*`a' + _b[dm2]*2.302585))) - (normprob(_b[/cut2] - (_b[IV]*`a' + _b[dm2]*2.302585))) in `b'
				replace pred4 = 1-(normprob(_b[/cut3] - (_b[IV]*`a' + _b[dm2]*2.302585))) in `b'	
				
				local a = `a' + 1				
				local b = `b' + 1			
				}		
					
					
			*Get cumultative predicted probabilities 	
				replace pred2 = pred1 + pred2
				replace pred3 = pred2 + pred3
				replace pred4 = pred3 + pred4
			
			*gen figure
				twoway (area pred1 xvalue,  lcolor(black)) (rarea pred1 pred2 xvalue,  lcolor(black)) (rarea pred2 pred3 xvalue,  lcolor(black)) (rarea pred3 pred4 xvalue,  lcolor(black)), legend(order(1 "Not at all important" 2 "Little importance" 3 "Fairly important" 4 "Very important")) name(tempg1, replace) graphregion(color(white)) ytitle("Cumulative probability")  ylab(0(0.2)1, nogrid) xtitle("Percentage Demand") 
